Efficient method for simulating quantum electron dynamics 
under the time dependent Kohn-Sham equation 
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A numerical scheme for solving the time-evolution of wave functions under the time dependent 
Kohn-Sham equation has been developed. Since the effective Hamiltonian depends on the wave 
functions, the wave functions and the effective Hamiltonian should evolve consistently with each 
other. For this purpose, a self-consistent loop is required at every time-step for solving the time- 
evolution numerically, which is computationally expensive. However, in this paper, we develop a 
different approach expressing a formal solution of the TD-KS equation, and prove that it is possible 
to solve the TD-KS equation efficiently and accurately by means of a simple numerical scheme 
without the use of any self-consistent loops. 
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I. INTRODUCTION 



Since the innovative work on the density functional 
theory (DFT) Q and the Kohn-Sham equation many 
kinds of static or adiabatic quantum electronic phenom- 
ena have been investigated based on first principles. As 
an extension of the DFT to non-adiabatic dynamical phe- 
nomena, the time-dependent density functional theory 
(TD-DFT) has been developed [§§]. By using the TD- 
DFT, some excitation phenomena have been analyzed 
more accurately than by using the DFT §. However, the 
formulation of the TD-DFT is too complicated to solve 
the wave functions numerically in order to see electron 
dynamics directly. So a considerable approximate for- 
mula called the TD-Kohn-Sham (TD-KS) equation has 
been applied for the numerical simulations ^|J^]. 

The difficulty in numerically solving the TD-KS equa- 
tion is the treatment of the density-dependent Hamilto- 
nian. The wave functions and the Hamiltonian should 
always be self-consistent with each other. A fourth order 
self-consistent iterative scheme was proposed by O. Sug- 
ino and Y. Miyamoto ||. However, the use of a SCF-loop 
at every time-step is computationally expensive. 

In this paper, we propose a new formalism for the nu- 
merical solution of the TD-KS equation. Based it on, 
we prove that a simple formula without SCF-loops can 
solve the TD-KS equation with sufficient accuracy. We 
find that computational techniques previously de- 

veloped by us for the one-electron TD-Schrodinger equa- 
tion in real space and real time are also useful for the 
TD-KS equation. 



II. CONVENTIONAL METHOD 



The TD-KS equation is a mean field approach used 
for describing the time-evolution of the electron density 
p via one-electron wave functions tp n under an effective 



Hamiltonian Tt 



Of 



H\p,t]i> n {t) ; H[p,t] 



A 



N 



+ V[p,t] , 



(1) 



v[ P ,t\ = v int [ P ] + v ext (t) , P (t) = liM*)l s 



71=1 

Here, V[p, t] is an effective potential which represents the 
internal mutual interactions Vint [p] and the external time- 
dependent potential V c ^ t (t)- Throughout this paper, we 
use the atomic unit h = 1, m = 1, e = 1 for equations 
and values. 

Due to the time-dependence of the Hamiltonian, the 
solution of the TD-KS equation can be formally ex- 
pressed in terms of a time-ordering exponential operator: 

Mt) = Texp[-i f dt'H\p,t']\ V>„(0) ■ (2) 

There are many numerical methods for computing 
Eq. (||) . The simplest method discretizes the elapsed time 
t into small time slices At, and approximates Eq. (0) as 



yj n (t + At) «exp -iAtH[p,t] i/> n (t) 



(3) 



and it is computed using the Runge-Kutta method, or by 
the split operator technique: 



4> n (t + At) ~ exp 



iAt A 
~~2 



exp 



exp 



1 

iAt A- 



Mt) ■ (4) 



However, this is not sufficiently accurate, because it ig- 
nores the time dependence of the Hamiltonian during the 
small time slice, while the splitting reduces accuracy to 
an even lower level. 

Another well-known computational method for Eq. (|^) 
uses a Hamiltonian in the middle of the steps, 

tp n {t + At) ~ exp \-iAtH[p,t+ ip n (t) . (5) 
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Eq. (||) is also computed by the split operator technique: 
ip n {t + At) 





■iAt A- 




exp 


~2~~2 _ 


exp 



Ai Trr At, 
T ^,t + T -] 



exp 



iAt A 



1pn(t) ■ (6) 



Here, V[p,t + At/2] is estimated from an interpolation 
between V [p, t] and V[p,t + At] . Therefore, they have to 
be solved by a self-consistent loop. This scheme is accu- 
rate enough; however, it is computationally expensive to 
perform the SCF-loop at every time-step. 

III. FORMULATION 

To avoid the use of a SCF-loop, we first express the 
time-evolution of wave functions using a Taylor develop- 
ment in exponential form as 



00 At k d k 
iP n {t + = E —^.^n(t) = exp 



k=0 



k\ dt k 



At 



d_ 

dt 



1pn{t) 



(7) 



We consider a quantity f{{ip}, {ip*},t) which depends 
on wave functions ip and time t explicitly. The time- 
derivative of this quantity is expanded by the chain rule, 



dl^d± Sf dr Sf_ 
dt dt ' SiJj dt ' Sip* 



df_ 

<9i cx 

Here, we have used the following notation, 

N 



d± SJ_ 
dt ' Sip 



E / dr 



dip m (r) df 



dt dip m (r) 



(8) 



(9) 



and d/dt cx means an explicate-time-derivative operator, 
which operates only explicitly-time-dependent quantities. 

By substituting the TD-KS equation (Q) into Eq. (|), 
the time-differential is generally expressed as 

i i = ( w ^)-^-(^,w-4 :r +i£. do) 

For example, it operates a wave function ip n as 

^ = (n[p,m- — -(n[p,tm -^7 + ^ 

= H[p,t]ip n , (11) 

because ip n does not depend on tp^ and t explicitly. 
Another example regards density p, 



m 

because p also does not depend on t explicitly. 



(12) 



By substituting Eq. @ into Eq. (0), we can formally 
write the solution without employing the time-ordering 
operator as 



r/) n (t + At) = exp ^ Un\p, #) • 4t 
1 L dip 



Mt) • (is) 



However, it does not describe the algorithm of computa- 
tions. To show the way of computation of Eq.(|l3|), we 
decompose the exponential operator as, 



tpn (t + At) ~ exp 
iAt 



At d 



exp ■ 



exp ■ 
At 



(A^) 



exp ■ 



2 <9t cx 
_5_ 
Sip 
_5_ 

Sip ' 



(V[p,t]iP) 
iAt r 



(V[p,t]iP) 
- (AiPT 



s 

Sip* 



5 

Sip* 



s 

Sip* 



exp 



At d 



2 <9i c 



Mt) • (14) 



Equation (|l4|) is correct up to the second-order of At. 

To clarify the meaning of the exponential operator 
which contains the Laplacian appearing in Eq. (flj), we 
expand it in a Taylor development as 



iAt 



exp ■ 



[(Aip) 



5 

4 L v_1 ~ y Sp 
{iAt) k - 



(Aipy 



Sip* 



Ipn 



E 

fc=0 



k\4 k 



Ipn ■ (15) 



The first-term (k = 1) of the series operates ip n as 

(Aip)~-(AiP)*~]ip n = Aip n . 

The second-term (k = 2) operates as 

2 

Ipn 



(16) 



^■ip-^-w 



AlP n 



(AiP) 



SAip n 
Sip 



Sip n 



A-jrj- ■ (AiP) = AAlp r , 



Generally, 



1p n = A k 1p n 



(17) 



(18) 



Thus, we obtain the following identity 
iAt 



exp ■ 



^■^p-^-w. 



tpn 

iAt 



exp 



-A 



ip n . (19) 
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Similarly, we expand the exponential operator which 
contains the effective potential appearing in Eq. ([li]) as 



At 
exp —r- 



k=0 



Ipn ■ 

(20) 



The first-term (k — 1) of the series operates ip n as 



(V[ P) t]i>) ■ — ~ (V[p, tW ■ ^\ i> n = V[p, 



(21) 



The second-term (k — 2) operates as 



(F[p,#)~ -(V[p,tW ~ 1>/V„ 



V[p,t]V[p,tty n + ((V[p,t}i,) 



sv\pA 

5ip* 



Sip 

Ipn 



V[p,t]V[p,t}ip n +((V[p,t]<p)-i> 

-((vipMr-^ 6 -^ 



* S V\p,t] 
Sp 



Ipn 



Ipn 



= V[p,t}V[p,t}<P n . (22) 
Thus, we obtain the following identity: 



exp ■ 



4>n 



exp 



^V[p,t] ip n . (23) 
i 



Substituting Eq. (HI),© into Eq. ©, we obtain, 
[At d 



ijj n (t + At) ~ exp 
-At 



exp 



'-v{ P ,t} 



exp 



2 dt ex 
iAt A 
~Y ~2 



exp 
exp 



-\At A 
-~2~ 2~ 
At d 



2 at 



Mt) • (24) 



As a result, we obtain the desired formula: 



ip n (t + At) ~ exp 
exp 



riAt A 

2~y 



^(Vint[p'}+V ext (t+^) 



exp 



iAt A 
1T1F 



^„(t) . (27) 



Here, V cxt (t + At/2) is the external force in the middle 
of the steps. Meanwhile, p' in VJ nt [//] is not the density 
in the middle of the steps, but it is the density after the 
preceding operation, namely 



p'(r) 



N 

£' 

n=l 



exp 



iAt A 
~Y~2 



ipn(r,t) 



(28) 



Therefore, the formula ( |27|) can be explicitly computed 
without employing any SCF loops. 

The present non-SCF formula ( p7j ) is quite similar with 
the conventional non-SCF formula (|J) and the conven- 
tional SCF formula (||). However, in this paper, we 
have derived the formula based on the strict solution (|l3| ) 
by considering the time-dependence of the Hamiltonian, 
while the conventional non-SCF formula did not con- 
sider the time-dependence. We can easily show that the 
present non-SCF formula is as accurate as the conven- 
tional SCF formula by associating p' with p(t + At) as, 

p> = J2 \Mt)\ 2 + iy(ftyi - ^nyC) t + G(At 2 ) 



n=l 
N 

Ei 

n=l 



2 

At 



At dp 
Pit + ^) + 0{At 2 ) 



dt 

-<3(Ai 2 ) 



(29) 



Therefore, both the non-SCF formula and the SCF for- 
mula are correct up to the second-order of At. 



By the way, Vj n t [p] does not depend on time explicitly, 
because the density p does not depend on time explicitly 
as shown in Eq. (jl2]). Meanwhile, V c ^t(t) does depend on 
time explicitly, 



gVjntM 

dt rx 







dt rx 



^0 



(25) 



Therefore, the exponential of the explicit-time- 
derivative operator appearing in Eq. ( p4| ) affects only the 
external time-dependent potential V c ^Jt) as 



exp 



[At d 



2 dt r 



V cxt (t) = V cxt (t + ^) (26) 



IV. COMPUTATIONAL TECHNIQUE 

Computational techniques previously developed by us 
for the one-electron TD-Schrodinger equation |s|JTo[| are 
also beneficial for formula (^7|). We discretize the wave 
functions in real space, and use the finite element method 
for spatial derivatives. The only difference in the scheme 
for the TD-KS equation and TD-Schrodinger equation is 
the exponential of the effective potential: 

<(r) =exp[— V int [p]\ Mr) • (30) 
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By this operation, the phase of the wave functions is al- 
tered at each point, but the density p(r) is not altered. 
Therefore, we take the value of Vj„t[jo](r) as a constant 
during the computation, which is calculated just before 
the computation. 

It is quite easy to improve the accuracy of formula (|27| ) 
to the fourth order. The fourth-order accurate formula 
is given by Suzuki's exponential product theory [0] as 



tp n (t + At) ~ S 2 (sAt; t + (l- s)At) 
S 2 {sAt: t+(l- 2s)At) S 2 ((l - 4s)Ai; t + 2sAt) 
S 2 {sAt-t + sAt) S 2 (sAt; t) tp n (t) 

Here, s and S 2 (At;t) are given as 

s = 1/(4- n) 



(31) 



(32) 



S 2 {At;t) =exp 



[iAt A 
2~~2~ 
At 

exp 



—v[p'A 



exp 



iAt A 
~~2 



(33) 



where, p' is the density after the preceding operations. 



V. EXAMPLE 

In this section, we perform a simple simulation to verify 
the efficiency and accuracy of the present method. The 
model system we use here is a one-dimensional isolated 
system in which two electrons interact by a delta-function 
interaction under an oscillating electric field. The two- 
body wave function ty(xi, x 2 ;t) in this system obeys the 
following TD-Schrodinger equation, 



i-Qj_^(xi,x 2 ;t) = 



1 d 2 
2dx 2 



1 d 2 



+ aS(xi — x 2 ) 



+ (xi + x 2 )E Q sin(cu t) Vp(xi, x 2 ; t) , (34) 

where a is the coupling constant of the interaction, and 
E Q is an external electric field to perturb this system. 

We suppose that ^{x\ ,x 2 ;t) is expressed by a common 
one-electron orbital wave function ip(x,t) as 

ty(xi,x 2 ;t) = ip(xi,t)ip(x 2 ,t) 

4= (x(T, rxMi, 02) - X(i, °i)x(t 02)) ■ (35) 

Thus, the TD-KS equation is derived exactly, 



• 9 ,. \ 



1 d 2 

20X 2 

p(x,t) = \^j(x,t)\ 2 . 



ap(x, t) + xE Q sin(w £) tp(x,t) , 

(36) 

We use the following parameters for computation: 



Size of the system 
Number of grid points 
Mutual interaction 
External force 
Frequency 
Small time slice 
Total time steps 



L = 8.0 
N p = 64 
a = 0.5 
E = 1/64 
u = 1/8 
At = 1/16 
N t = 256k 



First, we compute the lowest eigen state of this system 
using the time-independent Kohn-Sham equation: 



Eijj {x) 



1 d 2 

2dx 2 



ap(x) ip (x) 



(37) 



We use this state as the initial state. 

Second, we compute the time-evolution using Eq. (J27 
Third, by Fourier transforming the time-fluctuation of 
the polarization, we obtain the spectrum of the scattered 
light as shown in Figure [l]. 



0.30 




0.05 0.10 0.15 0.20 0.25 0.30 
Energy [a.u.] 

FIG. 1. Spectrum of the scattered light. A sharp peak 
found at 0.125 is corresponding to the Rayleigh scattering. A 
sharp peak found at 0.261 is corresponding to the emission 
from the first excited state to the ground state, this energy 
includes many-body and non-linear effects. 

The peak appearing in energy lj = 0.125 comes 
from the injected light. The peak appearing in energy 
u) = 0.261 is expected to be the excitation energy be- 
tween the first excited state and the ground state. 

We have calculated the excitation energy by certain 
other methods: Method (A) solves eigen states by the 
non-TD-KS equation (p7|), method (B) modifies the re- 
sult of (A) by using RPA, and method (C) diagonalizes 
the non-TD-Schrodinger equation. The results are listed 
below: 

Excitation energies calculated by some methods 



(A) 


non-TD-KS eq. 




= 0.199 


(B) 


non-TD-KS eq. with RPA 


WRPA 


= 0.255 


(C) 


non-TD-Schrodinger eq. 


w Sch 


= 0.260 




TD-KS eq. 


0J 


= 0.261 
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We found the peak obtained by the present method, 
i.e., the TD-KS equation, reproduces fairly accurately the 
excitation energy calculated by means of the exact diago- 
nalization of the non-TD-Schrodinger equation. Namely, 
by solving the TD-KS equation, dynamical phenomena 
can be described more accurately than using the RPA as 
far as the effective Hamiltonian is correct. 

Next, to evaluate the error of the method, we esti- 
mate the error of the density p(x, T) at a specified time 
T = 256[a.u.]. 

Error = / dx \p(x, T) - p cxact (x, T) | , (38) 
Jo 

here the exact value p exac t(x, T) is prepared in advance 
by performing the same simulation on an extremely small 
time slice At = 1/256 [a. u.]. 

Figure || shows the errors on some time slices obtained 
by three methods: the present non-SCF method (^), the 
conventional non-SCF method (|J), and the conventional 
SCF method (g). 



10" 1 - present 2nd order non-SCF — ^- 
conv. 1st order non-SCF • 
' conv. 2nd order SCF * 




4 8 16 32 64 128 

1/At 

FIG. 2. Errors in the density obtained by three meth- 
ods on some small time slices. The conventional non-SCF 
method is accurate up to the first order of At, while the 
present non-SCF method and the conventional SCF method 
are accurate up to the second order of At. In this test case, 
the error of the non-SCF method is almost as same as that of 
the SCF method. 

All methods are accurate enough in this result. How- 
ever, the conventional non-SCF method is stable only 
within a specific short time span: e.g. T — 512 [a.u.] 
for all At in this test. Meanwhile, the present non- 
SCF method and the conventional SCF method are sta- 
ble even in a long time span: e.g. T — 64M[a.u.], 
At = l/16[a.u.] in this test. Therefore, these methods 
are suitable for long time span simulations. 

We have also tested the simulation using the present 
fourth-order non-SCF method (|3l| ) and the fourth-order 
SCF method proposed in the literature |(|. Figure (j^) 
shows the errors. Both errors are much less than those 
of the second-order methods. 




10 -10 



4 8 16 32 64 128 

1/At 

FIG. 3. Errors in the density obtained by the fourth-order 
methods. Both errors are roughly proportional to At 4 , and 
they are much less than those of the second-order methods. 
In this test case, the error of the non-SCF method is almost 
as same as that of the SCF method. 



VI. CONCLUSION 



We have proved that simulation of the wave function 
under the TD-KS equation can be performed by a simple 
scheme and that there is no need for the use of SCF-loops 
to maintain the self-consistency of the effective Hamil- 
tonian. Our proposed non-SCF method is competitive 
in accuracy with the SCF method, and also it is supe- 
rior in computational efficiency. We are convinced that 
our method is helpful for investigating non-adiabatic and 
non-linear quantum electrons dynamics. 
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